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£S) ' Chemical reaction systems with a low to moderate number of molecules are typ- 



ically modeled as discrete jump Markov processes. These systems are oftentimes 
simulated with methods that produce statistically exact sample paths such as the 
Gillespie Algorithm or the Next Reaction Method. In this paper we make explicit 
use of the fact that the initiation times of the reactions can be represented as the 
firing times of independent, unit rate Poisson processes with internal times given by 
integrated propensity functions. Using this representation we derive a modified Next 
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Reaction Method and, in a way that achieves efficiency over existing approaches for 
exact simulation, extend it to systems with time dependent propensities as well as 
to systems with delays. 
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Due to advances in the knowledge of cellular systems, where there are low to moderate 
numbers of molecules of certain species, there has been a renewed interest in modeling 
chemical systems as discrete and stochastic as opposed to deterministic and continuous. 
Because of the intrinsic stochasticity at this level, understanding of a given system is gained 
through knowledge of the distribution of the state of the system at a given time. As it is 
typically impracticable to analytically solve for the distribution of the state of the system at a 
particular time for all but the simplest of examples, simulation methods have been developed 
that generate statistically exact sample paths so as to approximate the distribution. The 
two most widely used exact simulation methods are the original Gillespie Algorithm 5 ' 6 and 
the Next Reaction Method of Gibson and Bruck. 7 
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In this paper, we will explicitly represent the reaction times of discrete stochastic chemical 
systems as the firing times of independent, unit rate Poisson processes with internal times 
given by integrated propensity functions. Such a representation is not novel and is called 
a random time change representation in the mathematics literature. See, for example, 
Refs. However, using such a representation in an explicit attempt to develop new 

simulation methods has a number of benefits that have seemingly not been explored in the 
chemistry literature. First, the representation will naturally lead us to a modified version 
of the Next Reaction Method.- Second, the modified Next Reaction Method will be shown 
to be the natural choice for simulating systems with propensities that depend explicitly on 
time (such as systems with variable temperature or cell volume). Third, we will be able 
to easily extend our modified Next Reaction Method to systems that allow delays between 
the initiation and completion of reactions in a manner that achieves efficiency over existing 
methods. More precisely, in our modified Next Reaction Method for systems with delays 
no random numbers or computations will be wasted (such as happens in the method of 
Bratsun et al.— and Barrio et al.— ) (see Section I VII) . and there will be no need for the 
complicated machinery of the method developed by Cai^ (see Section IVII) in the handling 
of the stored delayed reactions. We note that the ideas we use to develop our modified Next 
Reaction Method are analogous to the theories of generalized semi-Markov processes^ 1 ^ 1 ^ 
and stochastic Petri nets^, and can also be extended to develop new accurate and efficient 
approximate tau-leaping methods. 18 

The outline of the paper is as follows. In Section [III we briefly present the original Gillespie 
Algorithm. In Section II III we introduce our representation of the reaction times as the firing 
times of independent, unit rate Poisson processes with internal time given by integrated 
propensity functions. In Section IIVI we rederive the Next Reaction Method and derive a 
modified Next Reaction Method using the representation detailed in Section IIIII In Section 
[V]we consider systems with propensities that depend explicitly on time and conclude that our 
modified Next Reaction Method is the preferable algorithm to use in such cases. In Section 
|VT]we consider systems in which there is a delay between the initiation and completion of 
some of the reactions and develop a new algorithm for simulating such systems that is an 
extension of our modified Next Reaction Method. 
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II. THE GILLESPIE ALGORITHM 



Consider a system consisting of N > 1 chemical species, {X±, . . . , X^}, undergoing M > 
1 chemical reactions, each of which is equipped with a propensity function (or intensity 
function in the mathematics literature), ai~(X). For the time being, assume that the time 
between the initiation and the completion of each reaction is negligible. To accurately 
simulate the time evolution of the number of each species, X(t) = {Xi(t), . . . , X N (t)} 6 N> , 
one needs to be able to calculate 1) how much time will pass before the next reaction takes 
place (i.e. initiates and completes) and 2) which reaction takes place at that future time. One 
can then simulate statistically exact sample paths for the system of interest. The following 
assumption, sometimes called the fundamental premise of chemical kinetics, is based upon 
physical principles and serves as the base assumption for simulation methods of chemically 
reacting systems:- 



where o(At)/At — > as At — > 0. Based upon the assumption (p^), the time until the next 
reaction, A, is exponentially distributed with parameter ao(X(t)) = ^2^=0 a k{X (t)) and the 
probability that the next reaction is the kth is a,k{X (t)) / ao(X (t)) . These observations form 
the foundation for the well known Gillespie Algorithm.-^ 

Algorithm 1. (Gillespie Algorithm) 

1. Initialize. Set the initial number of molecules of each species and set t = 0. 

2. Calculate the propensity function, a^, for each reaction. 

3. Set a = Ylk=i a k- 

4. Generate two independent uniform(0,l) random numbers T\ and r-i- 

5. Set A = l/aoln(l/rl) (equivalent to drawing an exponential random variable with 
parameter a ). 

6. Find jj, G [1, . . . , M] such that 



which is equivalent to choosing from reactions [1, . . . , M\ with the fcth reaction having 
probability ak/a . 



dk(X(t))At + o(At) = the probability that reaction k 



(1) 



takes place in a small time interval [t, t + At) 
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7. Set t = t + A and update the number of each molecular species according to reaction 
/'• 

8. Return to step 2 or quit. 

We point out that the Gillespie Algorithm uses two random numbers per step. The first 
is used to find when the next reaction occurs and the second is used to determine which 
reaction occurs at that time. In Section [IV] we will demonstrate how the Next Reaction 
Method generates exact sample paths while only needing one random number per step. 



III. REPRESENTATION USING POISSON PROCESSES 

We will explicitly represent the reaction times of chemical systems as the firing times of 
Poisson processes with internal times given by integrated propensity functions.- 1 *^ Using 
such a representation allows us to consider the system as a whole via a stochastic integral 
equation as opposed to solely considering how to calculate when the next reaction occurs 
and which reaction occurs at that time. The benefits of such a representation will stem from 
the fact that the randomness in the model is separated from the state of the system. 

Let v k , v' k G N" be the vectors representing the number of each species consumed and 
created in the /cth reaction, respectively. Then, if Rk{t) is the number of times that the kth 
reaction has taken place up to time t, the state of the system at time t is 

M 

X(t)=X(0) + J2Rk(Wk-Vk)- (2) 
fc=l 

However, based upon the assumption (CQ), R k (t) is a counting process with intensity a k (X(t)) 
such that Prob(i? fc (t + At) - R k (t) = 1) = a k (X(t))At for small At. Therefore, 

R k (t)=Y k y\ k (x( S ))d S y (3) 

where the Y k are independent, unit rate Poisson processes. Thus, X(t) can be represented 
as the solution to the following equation: 

X(t) = X(0) + ( [ a k (X(s))ds) K - Vk). (4) 

k=i ^° / 

Note that the state of the system, X(s), and hence each propensity function a k (X(s)), 

is constant between reaction times. In Section [V] we will consider systems in which the 



propensity functions are not constant between reactions, such as arise due to changes in 
temperature or cellular volume. 

We make two points that are crucial to an understanding of how different simulation 
methods arise from equation First, all of the randomness in the system is encapsulated 
in the Yj.'s and has therefore been separated from the state of the system. Thus, since the 
system only changes when one of the Y^'s change, the relevant question of each simulation 
algorithm is how to efficiently calculate the firing times of each and how to translate that 
information into reaction times for the chemical system. Second, there are actually M + 1 
relevant time frames in the system. The first time frame is the actual, or absolute time, t. 
However, each Poisson process Yj. brings its own time frame. More specifically, if we define 
Tjfc(t) = Jq ak(X(s))ds for each k, then it is relevant for us to consider Yj.(Tfc(t)). We will 
call Tfc(t) the "internal time" for reaction k. 

Definition 1. For each k < M, Tk(t) = J* * ak{X(s))ds is the internal time of the Poisson 
process Y k of equation (j4j). 

We will use the internal times of the system in an analogous manor to the use of "clocks" 
in the theory of generalized semi-Markov processes.— iJ&lS 

We now formulate the Gillespie Algorithm (Algorithm [TJ in terms of equation At 
time t, we know the state of the system, X(t), the propensity functions, a,k(X(t)), and the 
internal times, Tfc(t). Calculating 1) how much time will pass before the next reaction takes 
place and 2) which reaction takes place at that future time is equivalent to calculating 1) 
how much time passes before one of the Poisson processes, Y^, fires and 2) which Y& fires at 
that later time. Combining the previous statement with the fact that the intensities of the 
Poisson processes are a,k(X(t)) yields one step of the Gillespie algorithm. Use of the loss 
of memory property for Poisson processes (which negates knowledge of the internal times 
Tfc(t)) allows us to perform subsequent steps independently of previous steps. 

Note that in the Gillespie Algorithm the firing times of the individual processes Y& were 
calculated by first finding the time required until any of them fired, and then calculating 
which reaction fired at that future time. In the next section we show how the Next Reaction 
Method and our modified Next Reaction Method first calculates when each of the Yj. fires 
next, and then finds the specific reaction that fires by taking the minimum of such times. By 
not invoking the loss of memory property (and, hence, differentiating themselves from the 
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First Reaction Method-), the Next Reaction Method and modified Next Reaction Method 
make use of the internal times Tk(t) to nearly cut in half the number of random variables 
needed per simulation. 

IV. A MODIFIED NEXT REACTION METHOD 

We again consider the system (j3j). At time t we know the state of the system X — X(t), 
the propensity functions a k = a k (X(t)), and the internal times Tk = T/,(t). We also assume 
that we know Atk, the amount of absolute time that must pass in order for the kth reaction 
to fire assuming that stays constant over the interval [t, t + Atk). Therefore, T}~ = t + Atk 
is the time of the next firing of the kth reaction if no other reactions fire first. Note that if 
t — and this is the first step in the simulation of the system (and so Tk = 0), finding each 
Atk is equivalent to taking a draw from an exponential random variable with parameter a&. 
Because we know Atk, the internal time at which reaction k fires is given by Tk + akAtk- In 
order to simulate one step, we now note that the next reaction occurs after a time period of 
A = minfc{Atfc}, and the reaction that fires is the one for which the minimum is achieved, 
/i say. Therefore, we may update the system according to reaction /i, update the absolute 
time by adding A and update the internal times by adding a^A to Tk, for each k. 

For the moment we denote t = t + A and the updated propensity functions by The 
relevant question now is: for each k, what is the new absolute time of the firing of Yk, 
Tk, assuming no other reaction fires first? For reaction //, we must generate its next firing 
time from an exponential random variable with parameter a&. For k 7^ \x we note that, in 
general, the new absolute firing times will not be the same as the old because the propensity 
functions have changed. However, the internal time of the next firing of Yk has not changed 
and is still given by Tk(t) + akAtk- We also know that the updated internal time of Yk is 
given by Tk(t) = Tk(t) + Aa^. Therefore, the amount of internal time that must pass before 
the kth reaction fires is given as the difference 

(T k (t) + a k At k ) - (T k (t) + Aa k ) = a k (At k - A). 

Thus, the amount of absolute time that must pass before the kth reaction channel fires, At k , 
is given as the solution to a k Aik = ak(Atk — A), and so 

At k = ^(Atk- A), 
ak 
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Thus, we see that 

r k = ^{At k - A) + t = zr((t + At fc ) - (t + A)) + 1 = ^{r k - t) + t. 

dk &k 0-k 

We have therefore found the absolute times of the next firings of reactions k ^ fi without 
having to generate any new random numbers. Repeated application of the above ideas yields 
the Next Reaction Method.- 

Algorithm 2. (The Next Reaction Method) 

1. Initialize. Set the initial number of molecules of each species and set t = 0. 

2. Calculate the propensity function, ak, for each reaction. 

3. Generate M independent, uniform(0,l) random numbers r k . 

4. Set r k = l/a k \n(l/r k ). 

5. Set t = minfcjrfc} and let be the time where the minimum is realized. 

6. Update the number of each molecular species according to reaction /x. 

7. Recalculate the propensity functions for each reaction and denote by a k . 

8. For each k ^ //, set r k = (ak/~a~k)(T k —t)+t. 

9. For reaction \i, let r be uniform(0,l) and set r M = l/a M ln(l/r) + 1. 

10. For each k, set a k = a k . 

11. Return to step [5] or quit. 

Note that after the first timestep is taken in the Next Reaction Method, all subsequent 
timesteps only demand one random number to be generated. This is compared with two 
random numbers needed for each step of the original Gillespie Algorithm (Algorithm 1). 
We also note that the Next Reaction Method was originally developed with the notion of a 
dependency graph and a priority queue in order to increase computational efficiency (see Ref. 
E for full details). The dependency graph is used in order to only update the propensities 
that actually change during an iteration (and thereby cut down on unnecessary calculations) 
and the priority queue was used to quickly determine the minimum value in Step [5j We 
have omitted the details of these items as they are not necessary for an understanding of 
the algorithm itself. However, we point out that the use of a dependency graph in order 
to efficiently update the propensity functions is useful in any of the algorithms presented in 
this paper, and not just to the Next Reaction Method. 
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We now present an algorithm that is completely equivalent to Algorithm [21 but makes 
more explicit use of the internal times T k . In the following algorithm, we will denote by Pk 
the first firing time of Y k , in the time frame of Y k , that is strictly larger than T k . That is, 
Pk = min{s > Tk : Y k (s) > Y(T k )}. The main idea of the following algorithm is that by 
equation (j3J the value 

At k = {l/a k ){P k -T k ) (5) 

gives the amount of absolute time needed until the Poisson process Yk fires assuming that 
afc remains constant. Of course, a>k does remain constant until the next reaction takes place. 
Therefore, a minimum of the different Atk gives the time until the next reaction takes place. 
Thus, if we keep track of Pk and Tk explicitly, we can simulate the systems without the time 
conversions of step [8] of Algorithm [2j 

Algorithm 3. (Modified Next Reaction Method) 

1. Initialize. Set the initial number of molecules of each species. Set t — 0. For each k, 
set Pk = and Tk = 0. 

2. Calculate the propensity function, a& for each reaction. 

3. Generate M independent, uniform(0,l) random numbers r k . 

4. Set P k = ln(l/r fc ). 

5. Set At k = (Pk - Tk)/a k . 

6. Set A = minfcjAtfc} and let At^ be the time where the minimum is realized. 

7. Set t = t + A and update the number of each molecular species according to reaction 
//. 

8. For each k, set Tk = Tk + a^A. 

9. For reaction /i, let r be uniform(0,l) and set P M = P M + ln(l/r). 

10. Recalculate the propensity functions, a k . 

11. Return to step [5] or quit. 

We note that Algorithms [2] and [3] have the same simulation speeds on all systems that 
we have tested. This was expected as the two are equivalent. However, as will be shown in 
the next section, Algorithm [3] extends itself to systems with time dependent rate constants 
in a smooth way, whereas Algorithm [2] does not. We also point out another nice quality 
of Algorithms [2] and [31 Suppose that a system is governed by equation (j3J) except that the 
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Ya/s are no longer Poisson processes. That is, we suppose that the reactions do not have 
exponential waiting times, but have waiting times given by some other distribution. To 
modify Algorithms [2] and [3] to handle such a situation, one solely needs to change steps @] 
and M in each so that the waiting times are drawn from the correct distribution. 



V. TIME DEPENDENT PROPENSITY FUNCTIONS 

Due to changes in temperature and/or volume, the rate constants of a (bio) chemical 
system may change in time. Therefore, the propensity functions will no longer be constant 
between reactions. That is, dk(t) = a,k(X(t),t), and the full system is given by 

X(t) = X(0) +VYfe ( f a k (X(s),s)ds) (u' h - u k ), (6) 
k=i \ Jo / 

where the Y k are independent, unit rate Poisson processes. We consider how to simulate 

system using the Gillespie Algorithm, the Next Reaction Method, and our modified Next 

Reaction Method. 

The Gillespie Algorithm. At time t we know the state of the system, X(t), and, until 
the next reaction takes place, the propensity functions dk(X(t),s), for s > t. When the 
propensity functions depended only on the state of the system the Gillespie Algorithm 
calculated the time until the next reaction by considering the first firing time of M time- 
homogeneous Poisson processes. However, we now need to calculate the first firing time of 
M time-inhomogeneous Poisson processes. It is a simple exercise to show that the amount 
of time that must pass until the next reaction takes place, A, has distribution function 

/ M rt+A \ 

1 - exp I- J t ak(X(t), s)dsj . (7) 

Note that X(t) is constant in the above integrals because no reactions take place within the 
time interval [t,t + A). Using equation ([7]), A is found by first letting r be uniform(0, 1) 
and then solving the following equation: 

M rt+A 

/ a k (X(t),s)ds = \n(l/r). (8) 

In Appendix [A] we show that the reaction that fires at that time will be chosen according to 
the probabilities a,f e (X(t),t+ A)/a , where a = Ylh=i a k(X(t),t + A). Solving equation (jSj) 
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either analytically or numerically will be extremely difficult and time consuming in all but 
the simplest of cases. 

The Next Reaction Method. We begin by considering the first step of the Next Reaction 
Method. At time t = 0, we need to know the first firing times of independent, inhomogeneous 
Poisson processes. Therefore, we calculate the time that the kth reaction channel will fire 
(assuming no other reaction fires first) by solving for from: 



where is uniform(0, 1). Equation (jH]) can be solved either analytically or numerically. Say 
that reaction « is the first to fire and does so at time t. It is clear that to calculate the next 
firing time of reaction u we will need to generate another uniform(0, 1) random variable r M 
and solve 



What is less clear is how to reuse the information contained in tu for k 4 u. 

w 

Proceeding as in Ref. 13, denote by F n ^ a the distribution function for the nth firing of a 
reaction, where a is some parameter of the function. Gibson and Bruck prove the following: 

Theorem V.l (Gibson and Bruck's generation of next firing time^). Let r be a random 
number generated according to an arbitrary distribution with parameter a n and distribution 
function F antn . Suppose the current simulation time is t n , and the new parameter (after a 
step in the system in which this reaction did not fire) is a n+ \. Then the transformation 



generates a random variable from the correct (new) distribution. That is, t* has the correct 
distribution of the next firing time. 

Gibson and Bruck demonstrate use of the above theorem on a system whose volume is 
increasing linearly in time. In this specific case, it is possible to find closed form solutions 
of the distribution functions and their inverses. However, in general, calculating the distri- 
bution functions and their inverses may be a difficult (or impossible) task. For this reason 
Gibson and Bruck conclude "In general, it (the above method) may not be at all practicable 
and it may be easier to generate fresh random variables (according to the new distribution 
function)."- 




(9) 





(10) 
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In a situation in which Theorem IV. II is not practicable, the following steps must be taken 
to move one timestep beyond time t. First, generate uniform(0, 1) random variables, r k , and 
then solve 

a k (X(t),s)ds = \n(l/r k ) (11) 

for the time of the next firing of reaction k. We have been forced to use the loss of memory 
property of Poisson processes and generate new random variables. Thus, we are really now 
performing the First Reaction Method- 

The modified Next Reaction Method. As above, we begin by considering the first step 
of our modified Next Reaction Method. At time t = 0, we set T k = and P k = ln(l/rfe), 
where each r k is uniform(0, 1). To find the amount of time that must pass before the fcth 
reaction channel will fire if no others do first we solve for At k from: 

a k (X(0),s)ds = P k -T k = P k . (12) 



Again supposing that reaction [i fires first at time t, we update T k = J Q * a k (X(0), s)ds for 
each k. In order to calculate At M we must generate a new uniform(0, 1) random number, r M , 
set P M = P^ + hi{l/r^) and solve 

a^{X{t),s)ds = P fl -T tl . 



For ^ /j we still know that P k is the internal time of the next firing of reaction k, and 
so the amount of absolute time that must pass, At k , before the kth firing is given as the 
solution to 

/t+At k 
a k {X{t),s)ds = P k -T k . (13) 

Therefore, by keeping track of the internal times P k and T k we have been able to easily 
calculate the next firing of each reaction without having to generate another random number. 
We point out that using equation (fTBl to solve for the next firing time is no more difficult 
than using equation ffTTj) to find the next firing time in the Next Reaction Method, but 
equation fill I) demanded the generation of a random variable. We also point out that if 
there are closed form solutions to the above integral, such as the case of linearly increasing 
volume, then this method becomes very efficient. Further, even in this case of time dependent 
propensity functions, the modified Next Reaction Method easily lends itself to situations in 
which the waiting times between reactions are not exponential (only the generation of the 
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Pit's changes). We conclude that our modified Next Reaction Method will be preferable to 
either the Gillespie Algorithm or the Next Reaction Method on systems with propensity 
functions that depend explicitly on time. 

VI. SYSTEMS WITH DELAYS 

We now turn our attention to systems in which there are delays, r k > 0, between the 
initiation and completion of some, or all, of the reactions. We note that the definition of r k 
has therefore changed and is no longer the next reaction time of the Next Reaction Method. 
We partition the reactions into three sets, those with no delays, denoted ND, those that 
change the state of the system only upon completion, denoted CD, and those that change 
the state of the system at both initiation and completion, denoted ICD. The assumption 
(Cp) becomes the following for systems with delays: 

a k (X(t))At + o(At) = the probability that reaction k 

initiates in a small time interval [t, t + At), 
where o(At)/At — * as At — » 0. Thus, no matter whether a reaction is contained in ND, 
CD, or ICD, the number of initiations at absolute time t will be given by 

number of initiations of reaction k by time t = Y k (^J ak(X(s))ds\ , (15) 

where the Y\~ are independent, unit rate Poisson processes. 

Because the assumption (TTJJ, and hence equation ffl5|) . only pertains to the initiation 
times of reactions we must handle the completions separately. There are three different 
types of reactions, so there are three cases that need consideration. 

Case 1: If reaction k is in ND and initiates at time t, then the system is updated by losing 
the reactant species and gaining the product species at the time of initiation. 

Case 2: If reaction k is in CD and initiates at time t, then the system is updated only 
at the time of completion, t + r k , by losing the reactant species and gaining the product 
species. 

Case 3: If reaction k is in ICD and initiates at time t, then the system is updated by 
losing the reactant species at the time of initiation, t, and is updated by gaining the product 
species at the time of completion, t + r k . 
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The system can be written in the following integral form 



X(t) = X(0) +Y, Y *([ 

k&ND 



a k {X(s))ds (u' k - v k 



+ Zk ( I ak ( X ( s ~ T i)) ds ) i v 'k ~ v k) (16) 
+ J2 Wk (f a k (X(s-T k ))ds)v k - J2 Wk (f 



a k (X(s))ds ) v k 



kelCD v>/u ' kelCD 

where each a k (s) = for s < 0, and the Y k s, Z k s, and W k s are independent, unit rate 
Poisson processes. 

We note that there are more potential cases than those listed above. For example, the 
delay times, r k , may best be described as a random variable as opposed to being fixed or 
there could be multiple completion times for a single initiation (implying things happen in 
some order). For the sake of clarity we do not consider such systems in this paper but point 
out that it is a trivial exercise to extend the results of this section to such systems. 



A. Current Algorithms 

Based upon the discussion above, we see that simulation methods for systems with delays 
need to calculate when reactions initiate and store when they complete. However, because of 
the delayed reactions, the propensity functions can change between initiation times. Bratsun 
et al.— and Barrio et al.— used an algorithm for computing the initiation times that is exactly 
like the original Gillespie Algorithm except that if there is a stored delayed reaction set to 
finish within a computed timestep, then the computed timestep is discarded, and the system 
is updated to incorporate the stored delayed reaction. The algorithm then attempts another 
step starting at its new state. We will refer to this algorithm as the Rejection Method. 

Algorithm 4. (The Rejection Method) 

1. Initialize. Set the initial number of molecules of each species and set t = 0. 

2. Calculate the propensity function, a k , for each reaction. 

3. Set a = Y,k=i a k- 

4. Generate an independent uniform(0,l) random number, r\, and set A = 1/ao ln(l/r!). 
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5. If there is a delayed reaction set to finish in [t, t + A) 

(a) Discard A. 

(b) Update t to be the time of the next delayed reaction, /i. 

(c) Update x according to the stored reaction fi. 

(d) Return to step [2] or quit. 

6. Else 

(a) Generate an independent uniform(0,l) random number r 2 . 

(b) Find /i G [1, . . . , m] such that 

fl-l fl 

k=l k=l 

(c) If /i G ND, update the number of each molecular species according to reaction 
//. 

(d) If /i G CD, store the information that at time t + r^ the system must be updated 
according to reaction \i. 

(e) If fi G ICD, update the system according to the initiation of fi and store that at 
time t + r M the system must be updated according to the completion of reaction 

/'• 

(f) Sett = t + A 

(g) Return to step [2] or quit. 

At first observation the statistics of the sample paths computed by the above algorithm 
appear to be skewed because some of the timesteps are discarded in step However, 
because the initiation times are governed by Poisson processes via ffTol) . we may invoke the 
loss of memory property and conclude that the above method is statistically exact. 

The number of discarded A's will be approximately equal to the number of delayed 
reactions that initiate. This follows because, other than the stored completions at the time 
the script terminates, every delayed completion will cause one computed A to be discarded. 
Cai notes that the percentage of random numbers generated in step S] and discarded in 
step |Sa] can approach 50%.— Cai then develops an algorithm, called the Direct Method for 
systems with delays, in which no random variables are discarded. We present Cai's Direct 



Method below, however we refer the reader to Ref. 
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for full details. 



The principle of Cai's Direct Method is the same as that of the original Gillespie Algorithm 
and the Rejection Method above: use one random variable to calculate when the next 
reaction initiates and use another random variable to calculate which reaction occurs at 
that future time. However, Cai updates the state of the system and propensity functions 
due to stored delayed reactions during the search for the next initiation time. In this way 
he ensures that no random variables are discarded as in the Rejection Method. 

Suppose that at time t there are ongoing delayed reactions set to complete at times 
t + Ti,t + T 2 , . . . ,t + T d . Define T = and T rf+1 = oo. According to Cai's Direct Method, 
in order to calculate the time until the next reaction initiates, we first ask if the reaction 
takes place before t + Ti. If so, we may perform the step. If not, we must update the system 
according to the completion of the reaction due to complete at time t + Ti, update our 
propensity functions, and ask if the reaction takes place between t + T% and t + T%. In this 
manner we will eventually find when the next reaction initiates. Following the lead of Cai, 
we first present a method used for generating A.— 

Algorithm 5. (A generation for the Direct Method for systems with delays) 

1. Input the time t and ao = J2k a k- 

2. Generate an independent uniform(0,l) random number r±. 

3. If no ongoing delayed reactions, set A = 1/ao ln(l/ri). 

4. Else 

(a) Set i = 0, F = 0, and a t = clqT\. 

(b) While F < n 

i. Set F = 1 — exp(— a t ). 

ii. Set i = i + 1. 

iii. Calculate the propensity functions a^it + 7]) due to the finish of the delayed 
reaction at t + Ti, and calculate ao(t + Tj). 

iv. Set a t = a t + a (t + Ti)(T i+1 - T;). 

v. If i > 1 update the state vector x due to the finish of the delayed reaction at 
t + T^. 

(c) EndWhile 

5. Set i — i — 1. 
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6. Set A = Ti - (ln(l - n) + a t - a (t + Ti)(T m - T 4 )) /a (t + T t ). 



7. Endlf 



Because T 1; . . . , T rf are needed to perform the simulation, Cai introduces adx2 matrix, 
Tstruct, whose ith row contains Tj and the index ^ of the reaction due to complete at 
time t + Tj. During a simulation, if we find that A G [Tj,Tj + i), we delete rows 1 through 
% of Tstruct and set Tj = 7} — A for all of the other delay times. Also, rows are added 
to Tstruct when delayed reactions are initiated in such a way that we always maintain 
Tstruct(i, 1) < Tstruct(i + 1, 1). We present Cai's direct method below. 

Algorithm 6. (Direct Method for systems with delays) 

1. Initialize. Set the initial number of molecules of each species and set t — 0. Clear 
Tstruct. 

2. Calculate the propensity function, a,k, for each reaction. 

3. Set a = Ylk=i a k- 

4. Generate A via Algorithm [5j If A G [Ti,T i+ i) update Tstruct by deleting rows 1 
through i and update the other delay times as described in the above paragraph. 

5. Generate an independent uniform(0,l) random number r 2 . 

6. Find \i G [1, . . . , m] such that 



where the a^'s and ao are generated in stepHJ 

7. If fi G ND, update the number of each molecular species according to reaction \i. 

8. If fi G CD, update Tstruct by adding the row [r M , /i] so that Tstruct(i, 1) < 
Tstruct(i + 1, 1) still holds for all i. 

9. If G ICD, update the system according to the initiation of fi and update Tstruct 
by adding the row [r M , fj] so that T struct (i, 1) < Tstruct{i + 1, 1) still holds for all %. 

10. Set t = t + A. 

11. Return to step [2] or quit. 

We note that the Direct Method will use precisely one random number to find each 
initiation time. In this way the Direct Method is more efficient than the Rejection Method, 




k=l 



k=l 
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which discards a A (and therefore a random number) each time a delayed reaction completes. 
However, the extra machinery built into the Direct Method in order to find A will slow the 
algorithm as compared with the Rejection Method. Therefore, it is not immediately clear 
which method will actually be faster on a given system. 

B. The modified Next Reaction Method for systems with delays 

We now extend our modified Next Reaction Method to systems with delays. Recall 
that the central idea behind the modified Next Reaction Method is that knowledge of the 
internal time at which Y k fires next can be used to generate the absolute time of the next 
initiation of reaction k. The same idea works in the case of systems with delays because the 
initiations are still given by the firing times of independent Poisson processes via equation 
( I15p . Therefore, if T k is the current internal time of Y k , P k the first internal time after T k 
at which Y k fires, and the propensity function for the /cth reaction channel is given by a k , 
then the time until the next initiation of reaction k (assuming no other reactions initiate or 
complete) is still given by At k = {P k — T k )/a k . The only change to the algorithm will be 
in keeping track and storing the delayed completions. To each delayed reaction channel we 
therefore assign a vector, s k , that stores the completion times of that reaction in ascending 
order. Thus, the time until there is a change in the state of the system, be it an initiation 
or a completion, will be given by 

A = min{At fc , s k (l) - t}, 

where t is the current time of the system. These ideas form the heart of our Next Reaction 
Method for systems with delays: 

Algorithm 7. (Next Reaction Method for systems with delays) 

1. Initialize. Set the initial number of molecules of each species and set t — 0. For each 
k < M, set P k = and 7^ = 0, and for each delayed reaction channel set s k = [oo]. 

2. Calculate the propensity function, a k , for each reaction. 

3. Generate M independent, uniform(0,l) random numbers, r k , and set P k = ln(l/r k ). 

4. Set At k = (P k - T k )/a k . 

5. Set A = minfcjAtfc, s k (l) — t}. 
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6. Set t = t + A. 

7. If we chose the completion of the delayed reaction /i: 

• Update the system based upon the completion of the reaction //. 

• Delete the first row of s^. 

8. Elseif reaction \x initiated and \x G ND 

• Update the system according to reaction /i. 

9. Elseif reaction /i initiated and n G CD 

• Update s M by inserting t + into in the second to last position. 

10. Elseif reaction \i initiated and /i G ICD 

• Update the system based upon the initiation of reaction \x. 

• Update s M by inserting t + into in the second to last position. 

11. For each k, set = + a^A. 

12. If reaction \x initiated, let r be uniform(0,l) and set P M = P M + ln(l/r). 

13. Recalculate the propensity functions, a^. 

14. Return to step H] or quit. 

We note that after the first step, the Next Reaction Method for systems with delays 
only generates one random variable for each initiation as opposed to the two generated in 
the Direct Method. Further, Algorithm [7] performs the updates in a way that uses every 
random variable that is calculated yet does not have the complicated machinery necessary 
in the Direct Method. We should therefore expect that Algorithm [7] will need less time in 
the simulation of chemical reaction systems with delays then either the Rejection or Direct 
Method. We also note that similar to our modified Next Reaction Method, Algorithm [7J 
extends easily to systems with time dependent rate constants, and non-exponential waiting 
times between initiations. 

C. Numerical examples 

Example 1. Consider the following system consisting of two reaction channels: 

R x : X t +X 2 ^ X 3 R 2 : X 3 % 0. (17) 
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The reaction channel R\ belongs to ICD and R2 belongs to ND. Therefore, we update 
X\ — X\ — 1 and X2 = X2 — 1 at the moment of initiation of R\, but only update X3 = X3 + 1 
after a delay. Following Cai,^ we chose c x = 0.001, c 2 = 0.001, X^O) = 1000, X 2 (0) = 1000 
and X 3 (0) = 0. We let the delay of R\ be T\ — 0.1 and simulated this system from time t = 
until t = 1. These values were chosen so that the number of initiations that have delayed 
completions is approximately 100% of all initiations. Therefore, nearly 50% of all steps of 
the Rejection Method will discard a random variable, thereby maximizing its wastefulness. 

We performed 10 4 simulations using each of the Rejection, Direct, and Next Reaction 
Method for systems with delays. The Rejection Method of Barrio and Bratsun took 179.5 
CPU seconds, the Direct Method of Cai took 167.2 CPU seconds, and the Next Reaction 
Method took 82.8 CPU seconds. Therefore, the Rejection Method took 7.4% more time 
than the Direct Method and took 116.8% more time than our Next Reaction Method for 
systems with delays while the Direct Method took 101.9% more time than our Next Reaction 
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Method. We note that we have not reproduced the results stated in Ref. 
Method was found to be 23% more efficient than the Rejection Method. In fact, when the 
Direct and Rejection Methods are programed in such a way that the differences in the codes 
reflects the differences in the algorithms, one typically finds that the difference in simulation 
times does not differ substantially. Considering that for this example nearly half of all 
random numbers generated by the Rejection method in order to calculate A are discarded 
(which is a maximum in waste for the Rejection Method, see Ref. 113'), the fact that the 
Direct Method is not substantially more efficient than the Rejection Method points out that 
the time used by the steps in the Direct Method in order to calculate A is not negligible as 
compared to the time needed to generate random numbers. 

Because the Rejection Method becomes more wasteful as the number of rejected A's 
increases, we will test the three algorithms on a system in which we can easily control the 
percentage of A's that are discarded. 

Example 2. We consider a simple model of gene transcription whose non-delayed version 

The model consists of three species: gDNA (NN), messenger RNA 



can be found in Ref. 
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(mRNA), and the catalytic TProt. NN is assumed to be in such abundant quantities as to 
be constant, so the model is completely determined by the state of the species mRNA and 
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TProt. There are four reactions allowed in the model: 

?„ • Clh+TPrnt, 

(18) 



R x : NN klTProt ) mRNA R 3 : TProt 



P 2 : mPiVA ^0 R A : TProt ^ 0. 

We suppose that reaction one belongs to CD and has a delay of r = 5. It is sim- 
ple to show that the mean value of the state of the system has an equilibrium value of 



(mRNA, TProt) = ((^1^3) / (^2^4) , k^/k^), and the mean values of the propensities of the 
reactions have equilibrium values of 



\x = kxTProt = kx^ X 3 = k. 



'•3 



A 2 = k 2 mRNA = k x ^ A 4 = k 4 TProt = k 3 . 

Therefore, the expected percentage of the initiations that have delayed completions can be 
approximated by 7, which is given by 

fei 



- h = = t fc 4 , 1 o , 

Ai + A 2 + A 3 + A 4 2^1 + 2^ 2 fi + l' K ' 



For the Rejection Method, the number of discarded A's will be approximately the number 
of initiations of delayed reactions. Therefore the Rejection Method becomes more wasteful 
as the percentage of the total reaction initiations that have delayed completions increases, 
and so we may expect to see that as 7 increases the Direct Method will become relatively 
faster as compared to the Rejection method. To test this we set fc 2 = 1, k 3 = 15, and k^ = 1 
so that 7 = (l/2)ki/(kx + 1). kx now acts as a parameter that can be changed in order to see 
the effect 7 has on the relative speeds of the two algorithms. We note that the parameters 
were not chosen for their biological relevance, but instead were chosen for experimental ease. 

For a series of fci's we computed the CPU time needed for the Direct Method, Rejection 
Method, and Next Reaction Method for systems with delays to simulate the above system 
10 4 times from time to time 30. See Figure [TJ We see that as kx increases, the Rejection 
and Direct Methods remain relatively close in terms of efficiency with the Rejection Method 
being slightly more efficient for smaller kx and slightly less efficient for larger kx- However, 
the Next Reaction Method for systems with delays (Algorithm [7]), is significantly more 
efficient than both for all kx- 
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FIG. 1: The above plot compares the speeds of the Rejection Method, Direct Method, and Next 
Reaction Method for systems with delays as the percentage of timesteps that are rejected in the 
Rejection Method, as parameterized by k\, increases. For different values of k\, each method was 
used to simulate the system (fT8j) 10 4 times. The plot above gives the CPU time needed for each 
method as a function of k\. We see that the Rejection and Direct Methods are nearly equivalent 
while the Next Reaction Method for systems with delays is significantly more efficient than both 
for all k\. 

VII. CONCLUSION 

By explicitly representing the reaction times of discrete stochastic chemical systems with 
the firing times of independent, unit rate Poisson processes with internal times given by 
integrated propensity functions we have developed a modified Next Reaction Method. We 
extended our modified Next Reaction Method to systems with delays and demonstrated 
its computational efficiency on such systems over the Rejection Method of Bratsun et al. 
and Barrio et al., and the Direct Method of Cai. Considering that many models of natural 
cellular processes such as gene transcription and translation have delays between the initi- 
ation and completion of reactions, and that the Rejection method appears to be the most 
widely used method for simulating such systems, we feel that this extension will be useful. 
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Also, as is pointed out in the text, our modified Next Reaction Method can be easily ex- 
tended to systems with non-exponential waiting times between initiations and is preferable 
to both the Gillespie Algorithm and the original Next Reaction Method for systems with 
propensities that depend explicitly on time. We feel that having a single, efficient simulation 
method applicable to such a broad range of chemical systems will prove to be a beneficial 
contribution. 
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APPENDIX A: UNFINISHED CALCULATION 

In Section |V] we showed that if a system has propensity functions that depend explicitly 
on time, then the amount of absolute time, A, that must pass after time t before any reaction 
fires has distribution function 

/ M rt+A \ 

1 - exp I - ^ J a k (X(t), s)dsj . 

where r is uniform(0, 1). We will sketch the proof of why the reaction that fires at that time 
will be chosen according to the probabilities a k (X(t), t+A)/a , where a = J2 fc=1 a,k(X(t),t+ 
A). 

Let H{r) = J2l=i It +V a k(X(t), s)ds. For j < M, let At,- be the amount of time that must 
pass after time t before the jth reaction fires. Let F denote the random variable min{A£j}. 
Then, conditioning on the fact that F = A and using the independence of the underlying 
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(Al) 



Poisson processes we have 

P(At k < At jd -L k \F = A) = limP(At fe < At jti #.\F e [A, A + 5)) 

S— >o 

= Um PjAh < Atj^ k , F e [A, A + 5)) 

™ P(F e [A, A + 5)) 

= P{At k € [A, A + 6),Atjj^ k > A + 5) 
~ ™ exp(-if(A)) - exp(-#(A + 5)) 

P(At fc e [A, A + 6)) n j#fc P(Atj > A + g 

™ exp(-H(A))-exp(-H(A + 6)) 

It is a simple exercise to show that for any j < M 

P{Atj > s) = exp (- aj(X(t), s)ds^j . (A2) 

Combining equations (lAlj) and (IA2B with an application of L'Hopital's rule gives the desired 
result. 
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Figures 

Figure 1. 




Captions 

Caption for Figure 1. 

The above plot compares the speeds of the Rejection Method, Direct Method, and Next 
Reaction Method for systems with delays as the percentage of timesteps that are rejected 
in the Rejection Method, as parameterized by k\, increases. For different values of k\, each 
method was used to simulate the system (ITS]) 10 4 times. The plot above gives the CPU 
time needed for each method as a function of k\. We see that the Rejection and Direct 
Methods are nearly equivalent while the Next Reaction Method for systems with delays is 
significantly more efficient than both for all k\. 
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